if(1)
  N = 32;

  h = (1/N);
  % V = 1 + rand(N,N);
  V = ones(N,N);

  global Tree;
  global mvcount
  disp('Factorization');
  Tree = five_setup(V,N,h);
  disp('Extraction');
  DiagExact = five_extract(N,Tree);
  DiagExact = reshape( DiagExact, N*N, 1 );


end

if( 1 )

  disp('Generating H');
  nx = N;
  ny = N;
  H = sparse(nx*nx, nx*nx);
  for i1 = 1 : nx
    for j1 = 1 : ny
      pos = nx*(i1-1) + j1;
      pos1 = nx * (i1-1) + mod(j1,ny)+1;
      pos2 = nx * (i1-1) + mod(j1-2,ny)+1;
      pos3 = nx * (mod(i1,nx)) + j1;
      pos4 = nx * (mod(i1-2,nx)) + j1;
      H(pos, pos1) = -1;
      H(pos, pos2) = -1;
      H(pos, pos3) = -1;
      H(pos, pos4) = -1;
      H(pos, pos) = 4;
    end
  end
  H = H / h^2;
  for i1 = 1 : nx
    for j1 = 1 : ny
      pos = nx*(i1-1) + j1;
      H(pos, pos) = H(pos, pos) + V(j1,i1);
    end
  end

  G = inv(H);

end

global svdcase

if(1)
  svdcase = 5;
  nLevel = 3;
  mvcount = 0;


  disp('Constructing H matrix');
  [DiagBlock, SampleList, RKTree, nDiagBlock] = ConstructHMatrix( N, nLevel );
end

if(1)
  % Compare the L^2 error for each block on each level.
  MeshInd = reshape(1:N*N,N,N);
  Gerr = cell(nLevel,1);
  for iLevel = 1 : nLevel
    Gerr{iLevel} = zeros(RKTree{iLevel}.nRKBlock,1);
    for iRKBlock = 1: RKTree{iLevel}.nRKBlock
      Source = RKTree{iLevel}.RKSource(:,:,iRKBlock);  
      Target = RKTree{iLevel}.RKTarget(:,:,iRKBlock);  
      SourcePos = ...
	reshape(MeshInd(Source(1,1):Source(1,2), ...
	Source(2,1):Source(2,2)), [], 1);
      TargetPos = ...
	reshape(MeshInd(Target(1,1):Target(1,2), ...
	Target(2,1):Target(2,2)), [], 1);
      SourceInd = RKTree{iLevel}.RKSourceInd(iRKBlock);
      TargetInd = RKTree{iLevel}.RKTargetInd(iRKBlock);
      Gexact = G(TargetPos, SourcePos);
      Gappro = RKTree{iLevel}.RKBlockUniU{TargetInd} * ...
	RKTree{iLevel}.RKBlockUniM{iRKBlock} * ...
	transpose(RKTree{iLevel}.RKBlockUniU{SourceInd});
      Gerr{iLevel}(iRKBlock) = norm(Gexact - Gappro, 2);
    end % iRKBlock
  end % iLevel
end
